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Abstract. We have computed 3-D displace- 
ments induced by atmospheric pressure loading 
from 6-hourly surface pressure field from NCEP 
(National Center for Environmental Predictions) 
Reanalysis data for all Very Long Baseline In- 
terferometry) and SLR (Satellite Laser Ranging) 
stations. We have quantitatively estimated the 
error budget our time series of pressure loading 
and found that the errors are below 15%. We 
validated our loading series by comparing them 
with a dataset of 3.5 million VLBI observations 
for the period of 1980-2003. We have shown that 
the amount of power which is present in the load- 
ing time series, but does not present in the VLBI 
data is, on average, only 5%. We have also suc- 
ceeded, for the first time, to detect horizontal dis- 
placements caused by atmospheric loading. The 
correction of atmospheric loading in VLBI data 
allows a significant reduction of baseline repeata- 
bility, except for the annual component. 

1 Introduction 

Atmospheric pressure variations which can reach 
50 mbar induce deformation of the solid Earth at 
the level of several centimeters, which certainly 
should be taking into account in routine analy- 
sis of space geodesy observations. Atmospheric 
loading effects are computed by convolving the 
surface pressure field with Green's functions (see, 
for example, Farrell (1972)). 

van Dam and Herring (1994) and van 
Dam et al. (1994) have computed atmospheric 
loading using a former atmospheric model from 
the NMC (National Meteorological Center) with 
2.°5 by 2.°5 spatial resolution and 12 hour sam- 
pling and applied to VLBI (Very Long Base- 
line Interferometry) and GPS (Global Position- 
ing System) site position data. Analysis of reduc- 
tion of variance of baseline lengths showed that 
atmospheric loading signal was clearly present 



in VLBI and GPS datasets. However, only 62% 
of the power of the signal computed using their 
model was found in the VLBI data. This unex- 
plained discrepancy and very high computational 
cost of pressure loading calculation are the rea- 
sons why modeling atmospheric pressure load- 
ing computed using the global numerical weather 
model did not come into practice of routine data 
analysis. 

There are several reasons which motivated us 
to re-visit this topic: the accuracy of VLBI 
geodetic observations has increased during the 
last ten years which has improved our ability to 
detect tiny crustal motions. The NCEP/NCAR 
(National Center for Environmental Predic- 
tions / National Center for Atmospheric Re- 
search) Reanalysis project (Kalnay et al. (1996)) 
provides us a continuous and uniform dataset of 
surface pressure with a 6 hour temporal resolu- 
tion on a 2.°5 by 2.°5 grid for the period 1948- 
now. Because of the computer and network im- 
provements, it is now possible to compute on a 
operational basis the series of the atmospheric 
loading and apply them in a routine data analy- 
sis process, for example, for the Earth orientation 
service. 

In section 2, we describe our model of the at- 
mospheric pressure loading and present the er- 
ror budget of the loading time series. In sec- 
tion 3, we validate the time series of site dis- 
placements induced by pressure loading by re- 
analyzing a dataset of 3.5 million of VLBI obser- 
vations for the period 1980-2003. We have esti- 
mated global admittance factors for vertical and 
horizontal components. We also computed the 
reduction of variance of baseline lengths in order 
to compare our results with van Dam and Her- 
ring (1994). Discussion and concluding remarks 
are given in section 5. Brief outlines of the at- 
mospheric pressure loading service are given in 
section 6. 
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2 Computation of atmospheric pressure 
loading 

2.1 Characteristics of atmospheric pressure 
loading 

Displacements at the Earth's surface induced by 
surface pressure loading are computed by con- 
volving Green's functions (Farrell (1972)) with 
the surface pressure field from NCEP Reanalysis 
(Kalnay et al. (1996)). 

We model the ocean response to pressure forc- 
ing as the inverted barometer (IB). It has been 
shown that this model is adequate for describing 
sea height variations for periods typically longer 
than 15-20 days (see, for example, Wunsch and 
Stammer (1997)). However the ocean response 
to atmospheric pressure forcing significantly de- 
viates from the IB hypothesis at higher frequen- 
cies (Tierney et al. (2000)). We also assume 
that enclosed and semi-enclosed seas respond to 
atmospheric pressure as a " non-inverted barome- 
ter" , i.e. pressure variations are fully transmitted 
to sea floor. 

Figure and |2 shows the time series and their 
power spectrum for the period 2000-2003 at the 
Hartrao and Wettzell stations. They are repre- 
sentative of mid-latitude and equatorial inland 
stations. 




10 7 1<r" 10" 5 10-" 



0.50 - 




2000 2001 2002 2003 
year 



10 -s ^ 

10" , / 




10" 7 10" s 10" 5 10" a 



Figure 1: Vertical (left) and north (right) dis- 
placement induced by atmospheric loading for 
Hartrao station. 

All station displacements show significant 
narrow-band diurnal (Si), semidiurnal (S2) and 
annual (Sa) signals. Displacements for low- 
latitude stations are characterized by a strong 



wide-band annual and semi-annual signal and 
relatively weak amplitudes for periods below 
10 days (except the Si and S2 peaks). 
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Figure 2: Vertical (left) and east (right) dis- 
placement induced by atmospheric loading for 
Wettzell station. 

Mid latitude stations show the opposite be- 
havior; the mid-latitude atmospheric circulation 
is determined by the circulation of low and high 
pressure structures with typical period of 5- 
10 days. These periods are also the limit of the 
validity of the IB assumption. 

2.2 Error budget 

We would like to quantitatively asses the error 
budget of our atmospheric loading series. We 
identified four major sources of errors: 1) errors 
in the Green's functions, 2) errors in the surface 
pressure field, 3) errors in the land-sea mask and, 
4) mismodeling of the ocean response to pressure 
forcing. 

The Green's functions are computed for a 
spherically symmetric, non rotating, elastic and 
isotropic (SNREI) Earth model, adopting PREM 
(Dziewonski and Anderson (1981)) elastic pa- 
rameters. Therefore, we neglect the effects in- 
duced by Earth's anelasticity and ellipticity. The 
differences between our Green's functions and 
Green's functions for an anelastic Earth model 
(see, for example, Pagiatakis (1990)) are typ- 
ically below 1-2%. The effect induced by the 
Earth's ellipticity is of the order of magnitude of 
the Earth's flattening, i.e. 0.3%. 

Another source of possible errors is the errors 
in the surface pressure field from NCEP Reanaly- 
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sis. The first way to estimate them is to compare 
NCEP pressure held with ground pressure mea- 
surements. Velicogna et al. (2001) compared the 
differences between this atmospheric model and 
direct barometric observations. We would like to 
quantitatively asses the error budget, but since 
the spatial coverage of barometric measurements 
is not homogeneous, we choose another way to 
characterize possible errors, and we estimated 
the differences between two NCEP numerical 
weather models: NCEP Reanalysis and NCEP 
Operational Final Analysis. We have computed 
the station displacements with the NCEP Opera- 
tional surface pressure held with a spatial resolu- 
tion of l.°0. The RMS of the differences between 
the 3-D displacements computed with both pres- 
sure fields is, on average, about 10%. However, 
these differences are larger in mountainous areas; 
this is due to the fact that the 2° 5, or even the 
l.°0 spatial resolution of the NCEP Reanalysis 
and NCEP Operational datasets is not sufficient 
to model the topography in mountainous areas 
and, therefore, the surface pressure variations. 

Since the 2°5 resolution of the NCEP Reanaly- 
sis data is not sufficient to correctly represent the 
coastline, we choose a land-sea mask with a 0.°25 
resolution derived from FES99 (Lefevre et al. 
(2002)) ocean tidal model. In order to estimate 
the errors induced by the land-sea mask, we have 
computed 3-D site displacements with two land- 
sea masks with 0.°25 and 0.°50 resolution. The 
differences between the two time series do not 
exceed 5%. However the differences between the 
loading estimates with the 2° 5 and 0.°25 land- 
sea masks can reach 10% for the vertical compo- 
nents and for the 30% horizontal components, 
even for an inland station like Wettzell (Ger- 
many) 500 km far from the coasts. 

In order to estimate the errors induced by the 
mismodeling of the ocean response to pressure 
at high frequency, we compared ocean bottom 
pressure variations, as well as the induced load- 
ing effects, from two runs of the CLIO (Coupled 
Large-scale Ice Ocean) global circulation model 
(Goosse and Fichelet (1999)). The hrst run is 
forced by atmospheric pressure, surface winds 
and heat-fluxes (de Viron et al. (2002)), and 
the other one is forced only by winds and heat- 
fluxes, i.e. assuming an IB response. The differ- 
ences between these two runs are, therefore, in- 
terpreted as the departure of the ocean response 
to pressure forcing to the IB hypothesis. We have 
found that the mean vertical and horizontal dif- 



ferences are respectively below 10% and 20%. As 
expected these differences can be larger for island 
stations or stations close to the coast. 

Tabic 1: Global error budget of atmospheric 
loading estimates. 



Error source 


rms 


Green's functions 


2 % 


Surface pressure field 


10 % 


Land-sea mask 


5 % 


Ocean response 


10 % 


Total 


15 % 



Table ^ gives a summary of the global error 
budget. Combining all sources of errors, we can 
evaluate the uncertainty of the site displacements 
induced by atmospheric pressure loading at 15%. 

3 VLBI data processing 

3.1 VLBI data processing 

The processing of VLBI observations of time 
delays mostly follows the procedure described 
and Sovers et al (1998) and outlined in the 
IERS Conventions (McCarthy, 1996). Semi- 
empirical IERS96 nutation expansion was used. 
We modeled ocean tidal loading with GOT00 
(Ray (1999)) model for diurnal and semi-diurnal 
waves, NA099 (Matsumoto et al. (2000)) model 
for long-period zonal waves, and equilibrium tide 
for the polar tides and 18.6 year wave. NMF 
Mapping functions developed by Niell (1996) 
were used for modeling tropospheric path delay. 

3.2 Analysis of global admittance factor 

In solution Gl we computed admittance factors 
of the vertical, east and north components of 
the atmospheric pressure loading time series dis- 
placements into VLBI time delay, averaged over 
all sites. 

These admittance factors entered the estima- 
tion model the following way: 



a dr a dr 

t = > Al a fl7 A l2 a i2 - h 

orn dr i2 



i=l 



(1) 



where r is VLBI time delay, is the topocentric 
vector of the modeled site displacements due to 
atmospheric pressure loading for the jth station, 
r,j is the topocentric site coordinate vector. The 
index i runs through up, east and north com- 
ponents of the topocentric radius- vector. The 
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unknown vector of the admittance factors A was 
estimated in a single weighted LSQ fit together 
with site position and velocities, source coordi- 
nates, Earth orientation parameters (EOP) as 
well as 1.2 million nuisance parameters. 

It can be shown that these admittance fac- 
tors can be interpreted as correlation coefficients 
between the true (unknown) site displacements 
and our loading estimates within the validity of 
the assumption that the time series of the mod- 
eled site position displacements caused by atmo- 
spheric pressure loading and unmodeled contri- 
bution to time delay are not correlated. We refer 
the reader to Petrov and Boy (2003) for more 
details. 

We performed 4 different solutions. In solu- 
tion Gl we treated the three-dimensional vector 
of the admittance factor as a global parameter, 
which is common for all sites. Then we passed 
the time series of the site position variations 
through a narrow-band filter and extracted the 
annual component of the pressure loading sig- 
nal. In the solution G2 we estimated admittance 
factors for these time series. In the solution G3 
we estimated the admittance factor to the atmo- 
spheric pressure loading series with the annual 
component filtered out. Results are presented in 
Table [21 We also performed another solution in 
which we estimated admittance factors for each 
individual station. 

The estimates of the admittance factors of the 
time series of atmospheric pressure loading with 
removed annual component are closer to unity. 
Annual signals in site positions may be induced 
by other unmodeled signals such as hydrologic 
loading or thermal antenna height deformation 
(Nothnagel et al. (1995)) which are not mod- 
eled in this study. For example, van Dam et al. 
(2001) showed, that the contribution of continen- 
tal water loading can reach several cm at seasonal 
time scales. 

Contrary to wide-band signals, independent 
narrow-band signals are almost always corre- 
lated. Thus, condition of validity of this test is 
violated when we consider the annual signal. 

Table 2: Global Admittance Factors from VLBI 
solutions. 



Solution Up East North 

Gl 0.95 ± 0.02 1.16 ± 0.06 0.84 ± 0.07 

G2 0.46 ± 0.09 1.08 ± 0.26 -0.89 ± 0.26 

G3 0.98 ± 0.02 1.21 ± 0.06 1.02 ± 0.07 



Although the mean admittance factors are 
very close to unity, the admittance factors for 
individual stations are not always close to unity. 
We can often link anomalous admittance with 
the site positions. Stations located in mountain- 
ous areas or in the vicinity of the oceans are 
usually characterized by low admittance factors, 
even sometimes negative values. This can be ex- 
plained, as we showed in section 2.2, by errors in 
the atmospheric pressure field on mountainous 
areas and errors induced by the mismodeling of 
oceanic response to atmospheric pressure load- 
ing. 

3.3 Analysis of reduction of variance coeffi- 
cients 

In order to compare our results with results of 
van Dam and Herring (1994), we made another 
three solutions; in solution Bl, we did not apply 
the atmospheric loading contribution, but we did 
it in the solution B2. In solution B3 we have ap- 
plied the contribution of the atmospheric pres- 
sure loading time series with the annual compo- 
nent filtering out. 

We used the reduction of variance coefficient R 
of baseline lengths as a measure of agreement of 
the atmosphere pressure loading displacements 
series with observations: 



where Act 2 is the difference between the mean 
square of baseline length residuals before and af- 
ter adding the contribution due to station dis- 
placements caused by the atmospheric pressure 
loading, and a m is the variance of the atmo- 
spheric loading. 

A linear model was fitted in the series with dis- 
continuities at epochs of seismic events for sev- 
eral stations. The weighted root mean square of 
residual baseline lengths was computed for all 
baselines with more than 100 sessions for Bl, 
B2 and B3 solutions. 69 baselines fitted this 
criterion. The coefficients of reduction of vari- 
ance were computed using baseline length vari- 
ances. The mean coefficient of reduction of vari- 
ance of the B2 solution with respect to the ref- 
erence solution Bl is 0.86 ± 0.04 . The mean 
coefficient of reduction of variance of the B3 so- 
lution with respect to the reference solution Bl 
is 0.92 ±0.04 . 

Van Dam and Herring (1994) found a reduc- 
tion of variance of R — 0.62 and R = 0.76 with- 
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out the annual component. Our estimates are 
significantly closer to unity than theirs. We at- 
tribute these difference partly to improvements 
in modeling atmospheric loading, and partly to 
improvements of geodetic observations. Van 
Dam and Herring used an older atmospheric 
model provided by NMC (National Meteorologi- 
cal Center) with the same spatial resolution and 
a 12 hour temporal sampling compared to the 
6-hours for the NCEP Preanalysis model. The 
NCEP Reanalysis is also a uniform and continu- 
ous dataset, whereas the NMC pressure field had 
several discontinuities related to changes in the 
model. We used a land-sea mask with a higher 
resolution (0.°25) for modeling the IB response. 

4 Discussion and Conclusion 

A priori estimates of the errors of our time se- 
ries of the site displacements induced by atmo- 
spheric pressure loading are less than 15%. The 
estimates of the admittance factors from the so- 
lution Gl demonstrate that on average 95% of 
the power of the modeled pressure loading signal 
presents in the data well within the error bud- 
get. It allows us to make an important con- 
clusion that in average our model of the at- 
mospheric pressure loading quantitatively agrees 
with observations. The existence of the signifi- 
cant discrepancy between the model of the dis- 
placements caused by the atmospheric pressure 
loading and observations reported by van Dam 
and Herring (1994) has not been confirmed. Ex- 
cept for the annual component, applying the at- 
mospheric loading model results in a reduction of 
the power of the residual harmonic site position 
variations. For the first time we have detected 
horizontal component of atmospheric pressure 
loading in VLBI observations. 

Although we have an excellent agreement be- 
tween the model of atmospheric pressure loading 
and the observations in average, the estimates of 
the admittance factors significantly deviate from 
unity for some individual stations, suggesting de- 
ficiency of the model for these sites. These sta- 
tions are located either close to the coastline or in 
mountainous regions. In the first case, improve- 
ment of atmospheric loading modeling requires 
to model the high-frequency ocean response to 
pressure forcing which significantly deviates from 
the IB assumption. In the other case, the spatial 
resolution of NCEP Reanalysis data (2.°5) is too 
coarse to model the topography and therefore at- 



mospheric pressure variations. 

We evaluate the effects of modeling atmo- 
spheric loading on EOP estimation. We made 
two solutions, one with applying the contribu- 
tions of the 3-D atmospheric loading displace- 
ments, the other without. The RMS differences 
in polar motion and the UT1 angle between these 
two solutions are typically about 100 prad, which 
is 2-4 smaller than the current EOP uncertain- 
ties. 

5 Service of the atmospheric pressure 
loading displacements 

The tests described above allowed us to make 
a conclusion that the model of the atmospheric 
pressure loading has passed validation tests. We 
have computed continuous atmospheric pressure 
loading series for all VLBI and SLR sites start- 
ing from 1976.05.05. We found that when the at- 
mospheric pressure loading series are computed 
for two sites separated by 0.°05, the differences 
in site displacements never exceed the 1% level. 
Therefore, in the case if several VLBI or SLR sta- 
tions are located within 3 km from each other, 
the atmospheric pressure loading was computed 
only for one station and considered common for 
all stations located within this site. 

The first epoch of the atmospheric pressure 
loading time series is three days before the first 
observation used in data analysis and the last 
epoch is three days after the last observation of 
the site in the case if all stations at the site has 
ceased operations. Thus, these epochs are differ- 
ent for each site. 

The series of the atmospheric pressure loading 
for the sites, which are considered as active, i.e. 
continuing observations, are updated every day. 
The file with surface pressure from the NCEP 
Reanalysis for the last year is retrieved by ftp, 
split into monthly sections and stored. The at- 
mospheric pressure loading time series for active 
sites are augmented if the surface pressure files 
contain data for the epochs for which the pres- 
sure loading has not been computed. The NCEP 
Reanalysis numerical weather model are updated 
daily with the time lag 3-7 days. Our atmo- 
spheric pressure loading time series also updated 
with the time lag 3-7 days. 

Starting from December 2002 the atmospheric 
pressure loading contribution is incorporated in 
the model of VLBI data reduction in all so- 
lutions of the Goddard VLBI group, including 
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operational EOP solutions. The atmospheric 
pressure loading time series are available on the 
Web at http://gemini.gsfc.nasa.gov/aplo/ 
for all everybody without restrictions. 
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